Extended supersolid phase of frustrated hard-core bosons on a triangular lattice 



oo 

o 

O 
(N 

Oh: 

CD 



o 

CD 

e 

I 



O 

o 



o 

00 

o 



X 



Fa Wang, Frank Pollmann, Ashvin Vishwanath' 

' Department of Physics, University of California, Berkeley, CA94720 
(Dated: September 9, 2008) 

We study a model of hard-core bosons with frustrated nearest-neighbor hopping (t) and repulsion (V) on the 
triangular lattice. We argue for a supersolid ground state in the large repulsion (V ^ \t\) limit where a dimer 
representation applies, by constructing a unitary mapping to the well understood unfrustrated hopping case. This 
generalized 'Marshall sign rule' allows us to establish the precise nature of the supersolid order by utilizing a 
recently proposed dimer variational wavefunction, whose correlations can be efficiently calculated using the 
Grassman approach,. By continuity, a supersolid is predicted over the wide parameter range, V > —2t > 0. 
This also establishes a simple phase diagram for the triangular lattice spin 1/2 XXZ antiferromagnet. 

PACS numbers: 



SupersoHdity, where superfluid and crystalline orders co- 
exist, have fascinated physicsits since they were first theoreti- 
cally proposedl 1]. Recent experimental results in ^HelQl that 
are still under active debate have led to renewed interest. Ex- 
perimental developments on a different front, in the realiza- 
tion of optical lattices in ultracold atomic systems, motivated 
a search for a lattice supersolid. One of the more promising 
candidates is a model of strongly interaction hard-core bosons 
on a triangular lattice. The model Hamiltonian reads 



{if 



(ij) 



(1) 

where annihilates (creates) a hard-core boson on site i 

and rii — blb^ are density operators. The model is equivalent 
to the XXZ spin- 1/2 Hamiltonian on the triangular lattice 



H = 



2 



(2) 



where 3=*= — (l/2)(s^±is^), s^^^^^ are the Pauli matrices and 
are related to bosons by sf — {2ni — 1), = b'' , s~ = b, 
and Jz — V, J± — ~2t. The discussion below will be largely 
in terms of the bosons, although we will sometimes switch to 
the equivalent spin description, when that is more natural. The 
t > case coiTesponds to the unfrustrated hard-core boson 
model with repulsive nearest-neighbor interactions. For this 
case, a variety of studies including large scale quantum Monte 
Carlo simulations indicate a supersolid phase for all 

V/t > 8.9, stabilized by an 'order by disorder' mechanism. 
The solid order is of the three sublattice (++-) type, where two 
sublattices have the same boson density. 

The case of frustrated hopping (t < 0) suffers from a sign 
problem in the occupation number basis, and its ground state 
has been a subject of conjecture for the last three decades. 
The coiTesponding spin model is just the XXZ antifeiTomag- 
net, which, in the large limit was at the center of the RVB 
spin liquid proposal of Fazekas and Anderson . Later semi- 
classical and small cluster numerical studies suggested mag- 
netic order ItI], and general arguments which apply to the 
phase structure of bipartite dimer models, to which this model 



can be mapped in the large limit, also indicate the same 
resultfll^. However, the precise nature of ordering has not 
been conclusively established. In this letter, we show how this 
problem can be tackled, which is summarized briefly in boson 
language below. Due to frustration, the ground states in the 
V oo limit is extensively degenerate. Within this ground 
state manifold, we demonstrate that the frustrated problem 
with t < can be mapped, via a nontrivial unitary transforma- 
tion, onto the unfrustrated one with t > 0. Since the latter is 
weU understoodil|i|l[il, many properties of the frustrated 
case can be immediately derived. Such a generalized 'Mar- 
shall sign' was conjectured earlier based on state enumeration 
and numerics flEU- Here we construct the explicit transfor- 
mation which proves this conjecture, and moreover utilize it to 
deduce properties of the frustrated model. Our unitary trans- 
formation is diagonal in the occupation number basis, which, 
combined with our knowledge of the unfrustrated model, al- 
lows us to argue that a supersolid state is realized for < < 
as well. The precise details of the superfluid phase ordering 
requires further calculation. This is carried out using a vari- 
ational wave-function approach lfl2ll recently introduced for 
the unfrustrated t/V = 0+ limit, which captures the essential 
aspects of supersolid order very well and has good variational 
energy as compared to the quantum Monte Carlo results. Ap- 
plying the unitary transformation, we obtain a variation wave- 
function for the frustrated problem. Properties of this wave- 
function, in particular the phase correlations, are then calcu- 
lated. The state is found to be a supersolid and the resulting 
structure of the long range order (LRO) is shown in FIG. S^. 
Surprisingly, the superfluid amplitude vanishes on one of the 
sublattices and hence superfluidity lives exclusively on the 
honeycomb lattice formed by the remaining two sublattices, 
on which the amplitude alternates in sign. Contrary to naive 
expectations, the superfluid amplitude on these sites exceeds 
the maximum superfluid amplitude of the unfrustrated case. 

Finally, with this information in hand, we propose a phase 
diagram for the entire t/V > parameter range. Note the 
point t/V = 1/2 coiTesponds to the spin-isotropic triangu- 
lar antifeiTomagnet, where the 120° state is established. This 
can be smoothly connected to the large V supersolid state 
derived here as shown in FIG. [U Supersolid order would 
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FIG. 1 : Schematic phase diagram for both unfrustrated (t > 0) and 
frustrated hopping {t < 0) with repulsive interactions (V > 0). The 
three arrows are order parameters s — (b'' + b,ib — ib\2n — 1) on 
the three sub-lattices. For t/V < —1/2 or t/V > 0.1 there is only 
superfluid LRO (XY spin order). The thick line -1/2 < t/V < 
0.1 is the region of supersolid order. t/V — —1/2 is the SU(2) 
symmetric antiferromagnet. 

then naturally be preserved over the wide parameter range 
< —t < V/2, in contrast to the unfrustrated case, where it 
is only present for t < V/10. The frustrated triangular lattice 
boson model therefore appears to be an appealing candidate 
for the realization of the elusive supersolid phase - experi- 
mental prospects are discussed at the end. Note, this is also 
a phase diagram for the spin 1/2 XXZ magnet, and the regime 
of proposed RVB phase of Fazekas-Anderson |0] is actually a 
particular spin ordered state. 

Strong Repulsion Limit and Generalized Marshall Sign: 
In the limit of ^ we can restrict the Hilbert space 
to a manifold of states which correspond to classical Ising 
ground states of the triangular antiferromagnet |13]. Every 
such Ising configurations S can be represented by a close- 
packed dimer configuration C on the dual honeycomb lattice. 
This is a two-to-one mapping because of the Ising Z2 sym- 
metry (particle-hole symmetry in the boson language). The 
Hamiltonian projected into this degenerate subspace, in- 
troduces dynamics which splits the degeneracy. Note, to first 
order in degenerate perturbation theory, only the hopping term 
Hf ~ + H.c.) plays a role, leading to the 

double-hexagon resonance in FIG. |2](a) with amplitude — t. 
The problem of the large repulsion limit is therefore related to 
finding the ground state of a quantum dimer model with such 
dimer resonances. We have already noted that the t > case 
is tractable by Quantum Monte Carlo methods since there is 
no sign problem. However, the problem of interest here is the 
case t < 0. there is a unitary transformation which changes 
the sign of every matrix element of Ht, the problem can be 
mapped to unfrustrated case. This is generically not possible 
but, within the restricted Hilbert space, this indeed happens, 
and the required unitary transformation is the following. Con- 
sider the lattice in FIG.|2](c) with 1/4 special edges marked as 
thick and green. One can check by inspection that any double- 
hexagon resonance will change the number of covered special 
edges by ±2. Therefore, if we define a unitary transformation 
on the dimer basis 

\C') =U\C) =i^^'^^'>\C), (3) 

where i = and Ns{C) is the number of special (green) 

edges covered by a dimer in the dimer configuration C, the 




FIG. 2: (Color online) (a): two double-hexagon resonance configu- 
rations Cij and Cij = Cji. Red thick bars denote dimers. (b): Kaste- 
leyn orientation and edge weights of the honeycomb lattice. Thick 
blue edges have weight z, others have weight 1. The green dash-line 
rhombus encloses the enlarged unit cell. x,y are the principal axis. 
We use the six sites on a thick-edge hexagon as the basis, labeled 
as 1, . . . , 6 as shown in the right-bottom corner, (c): special edges 
(thick green on the honeycomb) for the unitary transformation relat- 
ing the unfrustrated and frustrated case. Thin solid green bonds on 
the triangular lattice are dual to the special edges. 

sign of the Hamiltonian will be changed. The unitary transfor- 
mation does not change the energy spectrum nor correlations 
that are diagonal in boson density. Hence, thermodynamics 
- that only depends on energy eigenvalues- is unchanged, for 
eg. transition temperature and nature of transitions. However, 
off-diagonal correlations are affected. We can therefore im- 
mediately conclude that the ground state has the same three 
sublattice density modulation as the supersolid phase in the 
unfrustrated model. Moreover, it also has a finite compress- 
ibility, identical to that in the unfrustrated problem, since this 
can also be expressed as a density-density correlation func- 
tion. The latter strongly suggests superfluid long range order 
(a 2D bosonic phase with finite compressibility at zero tem- 
perature), and taken all together this points towards supersolid 
order for t/V = 0^ as well. In order to directly establish off 
diagonal long range order, and obtain more detailed quanti- 
tative information, we turn to a variational wavefuction ap- 
proach. 

Variational Wavefunction We denote the two Ising states 
related to the dimer state C as S[C] and S[C] and consider the 
following kind of wavefunctions, 

\^) = ^</.(C)|C) - 5] 0(C) • {\S[C]) + \S[C])) /V2 (4) 
c c 

where (j>{C) is the (complex) amplitude. In the dimer repre- 
sentation, the projected Ht corresponds to the double-hexagon 
resonance in FIG.|2](a). Only those dimer configurations with 
'resonatable' double hexagons appear in the Hamiltonian ma- 
trix elements. We denote by c^, that a particular dimer cov- 
ering has a resonatable double hexagon at the pair of adjacent 
plaquettes i, j, where i is the plaquette with two dimers. Un- 
der resonance Cy Cij = cji. However, the rest of the dimer 
configuration with this pair of plaquettes removed dij remains 
unchanged. Hence, the entire dimer configuration may be de- 
noted as Cij + dij. Note, a single dimer configuration may 
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have many representations in this notation - one for each res- 
onatable hexagon pair. The variational energy E' = {'i'\Ht\'^) 
is 



+ c.c. 



(5) 



<ij> di- 



where c.c. is the complex conjugate. Before considering the 
frustrated case in detail, we briefly review the variational 
wavefunction for the unfrustrated casefl^]. There, t > 0, so 
the matrix elements of Ht are all non-positive in the dimer 
basis. Thus the Perron-Frobenius theorem applies and the 
ground state can be taken to be everywhere positive. Hence, 
we get a normalizable wavefunction if 0(C) = ^ P{C) with 
P{C) taken as the probability of the dimer configuration C. 
Equivalently, one can assign positive weights W{C) to each 
dimer configuration C, then the probability P{C) — W{C) /Z, 

where Z = Ec ^^(C)- 

The central assumption that leads to tractable wavefunc- 
tions is the following. We assign edge weights Wab to all hon- 
eycomb lattice edges ah, and write the weight of dimer cover- 
ing C as W{C) = ricovcrcd (afc) ^^b- Interpreting as a ficti- 
cious Gibbs weight, this corresponds to a problem of hardcore 
dimers in an external potential. Powerful Grassmann variable 
techniques have been developed for this problem, which will 
allow us to calculate properties of these wavefunctions. 

Plug the ansatz 0(C) — \/P{C) into and using the fact 

Pic - 4-d ■ ■ ) 

that the ratio „, is independent of the configuration 



E = ~t 



E 



^ P{c,j)P(c,j) 



(6) 



dij ) is the net probability of 



where P(cy) = Y.d,, ^(^r 
the local configuration. 

The dimer number operator is Uab = (1 + sf s|)/2 where 
ab is the honeycomb lattice edge dual to the triangular lat- 
tice edge ij. The probability P{cij), P{cij) are the expecta- 
tion values (ni2n34n56n78n9,io), (n23n45n67n89nio,i), re- 
spectively. This can be evaluated analytically by the Grass- 
mannian integral method 1 14]. In the Grassmannian formula- 
tion, the dimer partition function is represented as an integral 
over Grassmannian variables rja defined on the honeycomb 
lattice sites, Z = j ex^{Y,a,b''la^abTlb l'^)]\a'na = ^i\A], 
where Pf [A] is the Pfaffian of the Kasteleyn matrix A j 1 sll . 
and Aab = +Wab if the Kasteleyn orientation is from a 
to h, or = —Wab if otherwise (see FIG. |2] (b)). The prob- 
ability P{cij) is calculated as an expectation value in the 
Grassmannian theory, the rule is to replace Uab by AabVaVb, 

then we get P(Cy) = Wl2W34W56W78t«9,10 iUl^^iVa 

Thus the variational energy ^ can be written as E 

rlO \ 

-.iVa) 



<u'> 



nlO 
.=1 



with ^10,11 = Wio,i- 



The ten-point correlator of anticommuting 77 can be Wick- 
expanded into a Pfaffian of a 10 x 10 antisymmetric ma- 



trix. 



Ila=lVa) = Pf[(?7a'7b)] = \/ (ict[{T]aVb)] 



1 ... 10. The above formula can be further simplified to the 
determinant of a 5 x 5 matrix exploiting the bipartiteness of 

the honeycomb lattice: dlaii '7a) = I det[(?7Q776)]|, a = 
1,3,..., 9; b ~ 2,4, ...,10.. This is much more effi- 
cient than the brutal-force Wick expansion used by Sen et 
al. [12], which allows us to evaluate more complicated cor- 
relation functions later in this paper The two-point correlator 
iVaVb) — {A^^)ba can now be evaluated by a Fourier trans- 
formation since the Kasteleyn matrix A has 2D translational 
symmetry. For the chosen Kasteleyn orientation and basis 
shown in FIG. |2] (b)), in the thermodynamic limit, the two- 
point correlator of the site a in unit cell (0, 0) and the site b in 
unit cell {x, y) is 



{'>la,(Ofi)Vb,{x,y)) — 



2tt p2tt 



[A-\k)]i,ae''^''^''+''''y^ 
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where a,b = 1, . . . , 6, and A ^ (k) is the inverse of the 6x6 
anti-hermitian matrix A{k), 



Aik) 



03X3 R{k) 
-RHk) 03x3, 



withi?(/c) 



1 
z 



z 



where e-c ~ e'*^^ , ey — e^^^ . As is shown in Ref. fl2l for 
t > this variational wavefunction has two local minima 
at z sa 0.9258 with energy per site E = -0.13774i, and 
z K, 1.073 with energy per site E = — 0.13762i, correspond- 
ing to the two supersolid states, (H ) and (OH — ) of the 

triangular lattice boson model [d, 0, [sf]- For the frustrated 
case, the wavefunction is obtained by unitary transformation, 
I*') = U\^) '}2^(})'{C)\C),he.nctthtv3x\ai\onalvjawe- 
function is 0'(C) — i^^'^'^^ ^/ P{C). The variational energy 
E = (^'|iJt|^'') is of course the same as in the unfrustrated 
case. In order to understand the two variational wavefunctions 
better, we shall calculate two point correlation functions, as- 
suming for simplicity that the two points i, j are on the same 
horizontal line, and j is on the right. 

Diagonal Correlations: Consider first the density-density 
correlation function (s^Sj). Draw a line from i to j and it will 
cut through an set of honeycomb lattice edges < ab >. If the 
number of edges with no dimers cut by this line is even, then 
sf — +1 and otherwise — —1. In terms of the dimer num- 
ber operator nab the s^s^ becomes a non-local string operator. 



(^^4) = ( n - 1)) 

<ab> cut by ij 



(7) 



Expand the product we get 2l-'~*l terms(|j' — i\ is the distance 
between j and i measured by the triangular lattice constant), 
each of which is the type of correlation functions evaluated 
before. Because these operators are diagonal in the dimer ba- 



sis. 



I*') 



Off diagonal Correlations: The square of the off-diagonal 
long range order (ODLRO) parameter (blbj) is slightly more 
complicated. In the dimer basis it describes the simultane- 
ous resonances of two hexagons (if i and j are not neighbors). 
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FIG. 3: (Color online)Two possible simultaneous double-resonance 
needed for calculating (bj bj ): Ci,Cj ^ d, Cj , and d , Cj «-> d , Cj , 
with even(odd) number of no-dimer edges cut by the line i + x,j — 
a; (dash line). 



However for this process to happen sf and s| must be oppo- 
site. Label the two local resonating configurations on hexagon 
by Cj(j)and C,;(j), there are two possibilities of this si- 
multaneous double-resonance, shown in FIG.[3j with opposite 
conditions for the edges cut by the line i + x, j — x, where x 
is the horizontal triangular lattice vector The even(odd) re- 
quirement can be enforced using the dimer number operators 
as [1 ± J|(2nafc ~ l)]/2, where the product is over all edges 
< ab > cut by the line i + i, j — i (see FIG. [3]for an example). 
Consider t > case first, we have 



Wi2W34W5eW7sWQ,ioWii,i2 



I <^"-12'^34'^56"78"9,10"-ll,12 ' [1 + ]^(2nah - l)]/2 

"56"78"9,10"-11,12 • [1 - ]^(2nah - l)]/2^ | 
12 

a=l 

(8) 



n 



where the is the product of edge weights of the 

twelve(12) edges around hexagons i and j. Note, this sim- 
ple form arises because the string J|(2ria6 — 1) cancels out. 
We will see shortly that in the frustrated hopping case, this 
does not happen. 

If distance between i and j is large, the 12-point corre- 
lator (HaLi Va) can be factorized into two 6-point corre- 
lators (Ha^i '^a) • {Yi^^rVa}- And we have the relation 

^u;i2 W34 "'56^23^45 W61 1 (IIq^i Va) \ = So fiter- 

ally we have the factorization property 



^\b]\^')(^\b,m 



00 



For t < case we need to take care of the phases of cj)'. 
From FIG.|2](c) we can see that (^E*' | ^P') has similar form 
as the first line of (O, only the first term inside { • } acquires 
a minus sign. Therefore we get 



mblb^\^')=- 



^23^45^61^89^10,11^12,7 



U'12W34U'56W78'!«9,10^i'll,12 
?^12«34«56"78"9,10"■ll,12 J^lSriah ^ 1) 




n-l/2 = 



FIG. 4: Supersolid LRO from the variational wavefunction in the 
strong repulsion limit. Greyscale shows density order {2ni — 1) while 
arrows denote superfluid order {b]), for (a): Frustrated hopping {t < 
0) - note the sign structure of superfluid order; and (b):Unfrustrated 
hopping (t > 0) 



The product can be expanded into 2l^^'l^^ terms, each of 
which can be evaluated as before. Note, this correlation func- 
tion cannot be factorized as in the unfrustrated case and one 
necessarily needs to evaluate a string correlator. 

Results: We evaluate the above mentioned correlators up to 
distance \ j — i\ — 18 and extrapolate to infinite distance limit 
to determine the long range order. 

At the global energy minimum z = 0.9258 and 
for the unfrustrated (t > 0) case, the long range or- 
der parameter (s) = (fot + b,ib — ib\2n — 1) = 
(0.163,0,0.764), (0.372,0,-0.412), (0.372, 0, -0.412) for 
the three sublattices A,B,C, respectively (we have set the 
superfluid phase to zero and sublattice A is surrounded by 
weight z hexagon in FIG. |2] (b)). These numbers are in 
agreement with Quantum Monte Carlo (QMC) results. The 
average density deviation from 1/2 is |0.764 — 0.412 — 
0.4121/2/3 = 0.010, which is about 2%, in good agreement 



with QMC 1 4]. The solid order parameter is \nA + nBe 



2Tri/3^ 



jice 



47ri/3|2 



/9 = 0.0384 (nA,B,c boson densities on the 
three sublattices), which is about 15% smaller than the QMC 
result of 0.045 fl, [l^, but in good agreement with classical 
Monte Carlo calculations result 0.0389 of the same type of 
variational wavefunctions jll | . 

In the frustrated (t < 0) case the three sublattice order 
is (0,0,0.764), (0.389,0,-0.412), (-0.389, 0, -0.412), as 
shown in Fig. FIG. ID The average density deviates from 
1/2 by the same amount as the frustrated case. In spin lan- 
guage this means a non-zero average z-component of spin, 
I Sl/N\ = I Y.i sf/(2iV)l = 0.01, which is about 50% 
smaller than harmonic spin-wave result 0.02, which has the 
same symmetrylTJ. Note that surprisingly the superfluid 
amplitude(Xy-component of s) on the B,C sublattices is 
larger than those in the unfrustrated case, while it vanishes 
on the A sublattice. Note, this quantity can be directly mea- 
sured in Quantum Monte Carlo simulations of the unfrus- 
trated system in the large repulsion limit, by calculating cor- 
relations of the unitarily transformed operator. For example, 
with O = s^s~, with j to the right of i in the same hori- 
(9) zontal line and |j — i| > 2, the correlator to be measured is 



WOU = 



-X J 



Finally, we combine the present 
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results in the large repulsion (or V ^ — i) limit with known 
120° order in the isotropic V — ~2t ( — J±) limit. These 
can be connected without a phase transition, is we assume that 
the supersolid phase persists with no change in symmetry over 
the entire range V > —2t > 0. This scenario, which is also 
a phase diagram for the spin 1/2 XXZ antiferromagnet, is de- 
picted in Figure [T] (l^ 

Experimental Realization: How can the frustrated boson 
hoppings be experimentally realized? In lattice cold atom sys- 
tems, a recent experiment [16] demonstrated that 'repulsively' 
bound molecular bosons have frustrated hoppings. Consider 
preparing an initial state composed of molecules of pairs of 
atoms (either bosons or fermions) with one or zero molecules 
per site. If the interactions between atoms are now made re- 
pulsive, the effective molecular hopping is readily seen to be 
frustrated, since the singly occupied sites are lower in en- 
ergy. If this metastable state is sufficiently long lived, the 
equilibrium properties of the frustrated boson system can be 
accessed. In Josephson Junction Arrays, external magnetic 
fields can generate frustrated hopping lilTII . 
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